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ABSTRACT 

In  support  of  Task  07/040  (support  to  AIR  7000),  the  generation  of  correlated 
sea  clutter  returns  using  the  Memoryless  Nonlinear  Transform  is  investigated. 

The  generation  of  such  clutter  is  critical  to  the  performance  analysis  of  radar 
detection  schemes  under  realistic  clutter  scenarios.  This  feature  must  be  in¬ 
corporated  into  clutter  models  being  built  at  DSTO  to  test  radar  detector 
performance.  Examples  of  the  transform’s  application  is  given  for  a  number  of 
target  distributions  of  interest. 
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Simulation  of  Statistical  Distributions  using  the 
Memory  less  Nonlinear  Transform 

Executive  Summary 

This  work  supports  the  modelling  and  simulation  requirements  for  Task  07/040  (support  to 
AIR  7000).  Specifically,  in  order  to  test  and  evaluate  the  performance  of  radar  detection 
schemes  being  considered  at  DSTO,  it  is  necessary  to  evaluate  their  performance  in  a 
correlated  clutter  setting.  Models  of  interest  for  sea  clutter  correspond  to  those  being 
investigated  for  high  grazing  angle  detection,  specifically  the  K-  and  KK-Distribution 
models.  Realistic  sea  clutter  returns  exhibit  correlations  in  time,  and  it  is  the  purpose 
here  to  show  how  to  generate  correlated  K  and  KK-distributed  clutter  by  simulation. 

The  Memoryless  Nonlinear  Transform  (MNLT)  provides  a  general  framework  which  per¬ 
mits  the  generation  of  correlated  sea  clutter  returns  with  desired  marginal  distributions. 
Although  distributions  that  ht  into  the  class  of  spherically  invariant  random  processes 
(SIRV)  can  be  simulated  easily  without  the  MNLT,  some  distributions  of  interest  to  DSTO 
may  not  fit  within  the  class  of  SIRVs.  Hence  it  is  important  to  have  an  accurate  description 
of  the  MNLT. 

It  is  also  shown  how  shortcomings  of  the  MNLT  can  be  overcome.  Specifically,  it  has 
been  reported  that  the  MNLT  is  problematic  because  it  does  not  guarantee  an  exact  form 
of  desired  output  autocovaiance.  Essentially,  the  MNLT  converts  a  correlated  Gaussian 
process  into  a  new  process  with  desired  marginal  distributions.  It  is  also  correlated,  but  not 
through  the  same  autocovariance  function  as  the  Gaussian  process.  It  is  in  fact  related  to 
the  Gaussian’s  autocovariance  by  a  non-linear  mapping.  In  a  practical  situation  we  would 
like  to  specify  a  desired  process  with  given  autocovariance  function.  It  will  be  shown 
how  this  can  be  achieved,  using  an  inversion  method.  This  leads  to  a  novel  solution  to 
the  generation  of  correlated  returns,  with  desired  marginal  distributions  and  prescribed 
autocovariance  function. 

The  report  outlines  both  the  theoretical  aspects  of  this  transform,  and  how  it  is  imple¬ 
mented  in  practice.  A  simulator  has  been  developed,  which  is  fully  described  and  docu¬ 
mented  in  this  report.  A  number  of  simulations  are  included  to  show  how  the  code  can 
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be  utilised.  The  simulator,  written  in  Matlab,  is  in  a  form  that  will  enable  incorporation 
into  a  number  of  radar  models  under  development  at  DSTO. 
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1  Introduction 

The  simulation  of  sea  clutter  returns  is  an  integral  part  of  the  modelling,  analysis  and  val¬ 
idation  of  radar  systems.  In  particular,  it  is  important  to  have  the  capability  to  generate 
synthetic  radar  clutter  that  not  only  posesses  the  desired  marginal  distributional  charac¬ 
teristics,  but  also  exhibits  what  can  be  termed  short  range  correlations  in  time.  By  the 
latter  we  refer  to  clutter  returns  whose  covariances  decrease  with  time,  such  as  through, 
for  example,  a  damped  sinusoidal  function. 

The  purpose  of  this  report  is  to  show  how  the  Memoryless  Nonlinear  Transform  (MNLT) 
[1]  can  be  used  to  simulate  amplitude  distributions  with  correlation  in  time.  The  MNLT 
is  a  nonlinear  mapping  from  a  specified  probability  distribution  to  a  target,  or  desired 
distribution.  Its  formulation  means  we  have  the  flexibility  to  start  with  any  distribu¬ 
tion  and  convert  it  into  a  desired  distribution.  As  such,  we  can  begin  with  a  correlated 
Gaussian  stochastic  process,  and  produce  a  correlated  process  having  the  distribntional 
characteristics  of  interest. 

The  literature  contains  a  number  of  cases  where  the  MNLT  has  been  used  to  generate 
correlated  returns.  One  of  the  earliest  references  is  [2],  who  uses  the  MNLT  to  correlate 
Log-Normal  clutter.  Additionally,  [3]  also  examines  the  Log-Normal  case,  while  [4]  and 
[5]  extend  the  MNLT  to  the  Weibull  case.  The  extension  of  the  MNLT  to  the  simulation 
of  correlated  K-Distributed  clutter  first  appeared  in  [6],  followed  by  the  derivation  of  the 
relationship  between  correlation  fnnctions  in  [1]. 

In  the  DSTO  report  [7]  it  is  argued  that  the  MNLT  does  not  provide  an  entirely  sufficient 
solution  to  the  simnlation  of  correlated  clutter^.  Part  of  the  criticism  of  the  method 
is  that  it  does  not  indicate  how  the  correlation  characteristics  of  the  target  distribution 
can  be  controlled.  This  is  because  we  begin  with  a  correlated  Ganssian  process,  and 
produce  an  output  process  with  a  correlation  function  that  has  been  transformed  from  the 
Gaussian  correlation  function  via  the  MNLT.  However,  it  has  been  shown  in  [I]  that  the 
correlation  function  of  the  target  distribution  can  be  expanded  in  terms  of  the  Gaussian 

^It  is  worth  observing  that  [1]  appeared  after  the  DSTO  report  [7]  and  so  the  author  of  the  latter  was 
unaware  of  the  former. 
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correlation  function.  Specifically,  the  expansion  is  in  terms  of  Hermite  polynomials  from 
mathematical  physics.  This  relationship,  as  described  in  [1],  can  in  principle  be  used  to 
determine  which  Gaussian  correlation  should  be  used  to  produce  the  desired  correlations 
in  the  target  distribution  sample.  The  issue  with  the  relationship  is  that  it  expresses 
the  output  autocovariance  as  a  function  of  the  Gaussian  process  autocovariance.  What 
is  required  is  the  inverse  mapping,  which  will  then  specify  which  Gaussian  process  one 
should  start  with.  The  novelty  in  the  work  presented  here  is  that  a  general  methodology, 
supported  by  Matlab  code,  is  developed  to  achieve  this  aim.  It  involves  an  approximate 
inversion  of  the  correlation  mapping  from  [1],  which  then  enables  the  determination  of  an 
appropriate  Gaussian  input  process. 

In  addition  to  the  preceived  issues  with  determining  which  Gaussian  process  to  use  as 
input,  it  is  also  argued  in  [7]  that  the  theory  of  spatially  invariant  random  processes 
(SIRP)  provides  a  better  solution  to  the  problem  under  consideration.  While  this  is  true 
for  the  class  of  SIRPs,  the  problem  is  that  there  are  distributions  of  interest  to  DSTO 
that  may  not  fit  into  this  class.  The  prime  example  is  the  KK-Distribution  [20] ,  which  is 
a  model  for  high  grazing  angle  clutter.  At  the  time  of  writing  it  was  not  clear  whether 
the  KK-Distribution  fits  within  this  class.  Hence  it  is  important  to  examine  whether  the 
MNLT  is  useful  for  simulating  a  correlated  sequence  with  KK-Distributed  marginals. 

This  report  analyses  the  problem  of  generating  correlated  clutter  returns  from  a  mathe¬ 
matical  point  of  view,  but  also  develops  a  practical  mechanism  for  generating  it.  Included 
is  a  correlated  clutter  simulator  in  Appendix  B. 

To  clarify  what  is  the  purpose  of  this  work,  we  formulate  the  desired  outcome  in  the  fol¬ 
lowing  manner.  We  would  like  to  generate  a  random  sample  of  observations  xi,X2,  ■  ■  ■  ,Xn 
from  a  prescribed  continuous  random  variable  X,  with  density  function  fx{t),  cumulative 
distribution  function  Fxifi)  and  covariance  Co'v{Xi,Xj)  =  K{XiXj)  —K{Xi)K{Xj),  where 
the  random  variables  corresponding  to  the  realisations  Xi  and  Xj  are  Xi  and  Xj  respec¬ 
tively.  These  realisations  represent  amplitude  or  intensity  measurements.  We  refer  to  this 
stochastic  process  as  the  target  or  output  process,  since  it  is  the  result  of  the  MNLT.  The 
Gaussian  stochastic  process  used  to  generate  the  target  process  is  called  the  input  process. 
We  will  restrict  attention  to  stationary  processes  throughout,  so  that  the  mean  of  each 
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realisation  is  constant.  We  will  employ  the  notation  and  terminology  used  in  [1],  referring 
to  r{k)  =  K^XqX/:)  as  the  correlation  function  of  the  process.  We  furthermore  define  the 
covariance  coefficient  of  the  process  as  (r-(fc)  —  IE(Alo)E(lAfc))/Var(X). 

Such  a  series  of  observations  can  be  used  as  range  samples  from  a  single-pulse  return  that 
are  input  to  a  detection  scheme,  such  as  a  Constant  False  Alarm  Rate  (CFAR)  detector,  or 
as  time  samples  from  multiple  pulses  in  a  single  range  cell.  We  do  not  examine  long-term 
correlations  in  this  work. 

This  report  is  structured  as  follows.  Section  2  introduces  the  MNLT,  and  Section  3  ex¬ 
amines  the  relationship  between  correlation  functions,  and  shows  how  they  can  be  related 
analytically  in  a  special  case.  Section  4  discusses  our  implementation  of  the  algorithms 
in  question,  while  Section  5  discusses  some  examples  of  its  output.  Appendix  A  describes 
the  use  of  the  simulator  from  the  point  of  view  of  the  user,  while  Appendix  B  contains 
the  source  code  for  the  simulator.  Finally,  Appendix  C  contains  numerically  determined 
values  related  to  the  function  described  in  Section  3. 

Due  to  the  fact  that  this  work  is  highly  probabilistic  in  nature,  it  is  worthwhile  listing 
some  useful  references  on  statistics,  probability  and  simulation.  Suitable  references  on 
probability  and  stochastic  processes  include  [8]  ,  [9]  and  [10],  while  [11]  is  a  useful  reference 
on  simulation.  Statistical  tests  and  methods  are  described  in  [12].  Standard  radar  texts 
all  describe  sea  clutter  and  their  characteristics,  for  example  [13]  and  [14]  both  describe 
sea  clutter  in  depth. 
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2  Simulation  via  The  Memoryless  Nonlinear 

Transform 

The  MNLT  is  essentially  an  extension  of  a  fundamental  method  with  which  a  single  random 
variable  return  can  be  generated.  We  thus  begin  with  the  generation  of  independent 
returns,  and  show  how  this  can  be  modihed  to  produce  dependency  between  the  marginal 
distributions  in  the  process  being  simulated.  Throughout,  we  are  focusing  on  continuous 
random  variables  as  specified  in  Section  1.  We  also  employ  the  notation  X  =  Y  to  mean 
that  the  random  variables  X  and  Y,  defined  on  a  common  probability  space  (0,.^,  P), 
share  the  same  distribution  function.  This  means  that  VA  £  Fx{A)  =  P(X  €  A)  = 
Fy{A)  =  P(y  e  A). 

The  following  Lemma  is  the  hrst  key  simulation  result: 

Lemma  2.1  Suppose  X  is  a  continuous  random  variable  with  cumulative  distribution 
function  Fx,  and  let  R  =  R{0, 1)  be  a  uniformly  distributed  random  variable  on  the  unit 
interval.  Then  the  random  variable  Ff^^{R)  =  X. 

To  see  this,  let  Y  =  Ff^^{R)  and  observe  that  the  cumulative  distribution  function  of  T  is 

Fviy)  =  nFf^\R)  <y)=  ¥{R  <  Fx{y))  =  Fxiy),  (I) 

using  the  cumulative  distribution  function  of  a  uniform  random  variable.  This  implies 
Y  =  X,  as  required. 

Lemma  2.1  states  that  to  generate  a  realisation  from  a  distribution  X  we  need  only  generate 
a  uniform  random  number  between  0  and  1  and  evaluate  the  inverse  of  the  cumulative 
distribution  function  at  this  point.  Since  we  are  dealing  with  continuous  distributions,  the 
cumulative  distribution  function  will  be  monotonically  increasing  and  continuous,  and  so 
the  inverse  will  always  exist. 

To  illustrate  how  (1)  can  be  used  to  simulate  a  random  variable,  consider  the  case  where 
X  is  an  exponential  distribution  with  density  fx{x)  =  e~^,  for  a;  >  0.  Then  it  is  not 
difficult  to  show  that  its  cumulative  distribution  function  Fx{x)  =  1  —  e~^  and  Ff^^{x)  = 
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—  log(l  —  x).  Hence  Lemma  2.1  implies  that  if  r  G  [0,1]  is  a  random  number,  then 

—  log(l— r)  is  a  realisation  of  X.  Generating  a  large  sequence  of  such  numbers,  and  plotting 
a  histogram,  will  show  that  they  have  the  same  statistical  characteristics  as  realisations 
from  an  exponential  distribution.  It  is  worth  noting  that  if  r  is  a  random  number  between 
0  and  1,  so  is  1  —  r  and  so  for  simulation  purposes  we  can  focus  on  generating  —  log(r)  in 
the  above.  This  result  is  clarified  in  the  following  Lemma: 

Lemma  2.2  If  R  =  R{0, 1),  then  1  —  R  =  R{0, 1). 

Lemma  (2.2)  is  proven  easily  by  constructing  the  distribution  function  of  1  —  i?.  Another 
useful  result  is  presented  in  Lemma  2.3,  which  is  essentially  equivalent  to  Lemma  2.1: 

Lemma  2.3  If  X  is  a  continuous  random  variable  with  cumulative  distribution  function 
Fx  and  R  =  R{0, 1),  then  Fx{X)  =  R.  That  is,  the  random  variable  resulting  by  substi¬ 
tuting  X  into  its  distribution  function  is  uniformly  distributed. 

This  Lemma  is  proven  by  showing  that  the  distribution  function  of  Fx{X)  matches  that 
of  R. 

Lemma  2.1  can  be  used  to  generate  independent  samples.  Recalling  the  problem  specih- 
cation  from  Section  1,  we  want  to  generate  a  correlated  sequence,  whose  point  to  point 
distributions  match  a  prescribed  distribution.  It  is  clear  the  only  way  to  proceed  would 
be  to  generate  a  correlated  sequence  of  uniform  random  numbers,  and  then  apply  Lemma 
2.1  to  produce  the  sequence  with  target  marginal  distribution.  This  sequence  will  be  thus 
correlated.  However,  we  can  use  Lemma  2.3  together  with  Lemma  2.1  to  do  this  more 
generally.  The  key  to  this  is  the  following  version  of  Lemma  2.1: 

Lemma  2.4  Suppose  Xi  and  X2  are  two  random  variables,  with  distribution  functions 
Fx^  and  Fx2  respectively.  Then  the  random  variable  Ff^^{Fx^{Xi))  =  X2. 

This  Lemma  follows  by  applying  Lemma  2.3  to  Lemma  2.1. 

Lemma  2.2  is  the  key  to  generating  correlated  realisations  of  a  random  process.  Define 
a  function  ri{x)  =  Ff^^(Fxi{x)).  Then  Lemma  2.2  implies  that  if  x  is  a  realisation  of 
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iV(0,l) 


R{0,1) 


^  Target  Distribution 


Figure  1:  The  Memoryless  Nonlinear  Transform  converts  a  Gaussian  variable  into  one 
of  the  target  distribution. 


random  variable  Xi,  then  r]{x)  is  a  realisation  of  random  variable  X2.  Suppose  we  have  a 
wide  sense  stationary  stochastic  process  f{t)  evolving  over  time  t,  with  correlation  function 
r^{T)  =  ]E((^(0)C(t)),  and  each  (^{t)  =  Xi.  Then  we  can  transform  this  process  to  produce  a 
new  stochastic  process  9{t)  =  r]{C{t)),  which  will  be  distributed  according  to  X2  pointwise 
{0{t)  =  X2,  Vt  >  0).  However,  this  new  process  will  have  a  correlation  function  r0(r)  = 
K(r^(r)),  where  k(-)  is  a  nonlinear  function. 

The  utility  of  this  and  Lemma  2.2  is  that  we  can  generate  any  random  variable  X2  from  any 
Xi.  Hence,  we  choose  Xi  to  be  a  Gaussian  random  variable,  and  generate  a  correlated 
Gaussian  sequence.  This  is  a  relatively  straightforward  exercise,  see  for  example  [15]. 
This  is  then  transformed  into  a  sequence  of  random  variables  with  the  desired  statistical 
distribution  using  the  function  rj.  Figure  1  illustrates  the  MNLT  process  graphically. 

The  resulting  sequence  will  consist  of  dependent  realisations.  Due  to  the  nonlinear  nature 
of  77,  the  output  sequence  will  not  necessarily  have  the  same  correlation  properties  as  the 
Gaussian  input  process.  This  is  examined  in  the  next  Section. 
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3  Correlation  Function  Mapping 

This  Section  is  concerned  with  the  autocorrelation  properties  of  the  MNLT,  and  in  par¬ 
ticular,  investigates  the  relationship  between  the  Gaussian  input  and  target  distribution 
output  autocorrelation.  As  remarked  previously,  [1]  derives  an  expansion  of  the  output 
autocorrelation  in  terms  of  the  input  autocorrelation.  What  is  desired  is  an  inverse  map¬ 
ping,  so  that  we  know  which  Gaussian  process  to  begin  with,  in  order  to  generate  a 
desired  stochastic  process  with  specified  marginal  distributions  and  autocorrelation  func¬ 
tion.  While  this  is  not  always  possible  to  do  analytically,  it  can  be  done  through  a  series 
of  functional  approximations.  We  begin  with  the  series  representation  of  the  output  au¬ 
tocorrelation. 


3.1  Correlation  Function  Expansion 

In  order  to  understand  the  relationship  between  the  correlation  functions  related  to  the 
MNLT  it  is  necessary  to  take  a  brief  digression  into  some  special  functions  from  math¬ 
ematical  analysis.  A  useful  reference  on  the  latter  is  [16],  although  we  merely  require 
the  properties  of  certain  special  functions  for  the  work  to  follow.  Hermite  polynomials 
[16,  17]  are  a  useful  set  of  functions  used  in  mathematical  analysis,  and  in  particular, 
in  Hilbert  space  theory.  They  form  a  complete  orthogonal  sequence  with  respect  to  a 
Gaussian  weight  function,  and  arise  in  the  theory  of  differential  equations — in  particular, 
as  a  special  case  of  the  Sturm-Liouville  boundary  value  problem  [17,  18].  The  0th  order 
Hermite  polynomial  is  Ho(t)  =  1  and  for  n  £  N,  the  nth  order  Hermite  polynomial  is 

9  9 

=  (-l)"e*  — (e-^  ). 

Hermite  polynomials  can  also  be  generated  with  a  second-order  recurrence  relation.  Noting 
that  Hi[t)  =  2t,  it  can  be  shown  that 

Hn+i{t)  =  2tHn{t)  -  2nHn-i{t). 
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This  enables  us  to  sequentially  generate  Hermite  polynomials.  The  generating  function  of 
the  Hermite  polynomials  is 

OO  . 

n\ 

n=0 

It  is  worth  noting  that  the  generator  (2)  essentially  specifies  a  Gaussian  density  in  terms 
of  the  Hermite  polynomials. 

The  following  Lemma  is  an  interesting  result  that  will  be  applied  in  the  Lognormal  example 
considered  in  the  next  Section: 

Lemma  3.1  Suppose  X  =  A^(/r,0.5)  is  a  Gaussian  random  variable  with  fixed  mean  pL. 
Then  for  each  n  G  N,  E  {Hn{X))  =  (2/x)”. 


The  proof  follows  directly  from  the  generator  function  (2). 


We  now  derive  the  relationship  between  Gaussian  input  and  target  output  autocorrelation 
functions.  It  essentially  follows  from  the  generating  function  (2).  Using  equations  (9)  and 
(10)  of  [19],  (with  the  choice  of  j  =  1,  A:  =  1,  z  =  |,  x  =  -^  and  y  =  ^),  it  can  be  shown 
that 


e 


a/i  -  r^e 


(fi+yfi 


n=0 


(3) 


As  pointed  out  in  [1],  this  equation  gives  an  expansion  of  the  generator  of  an  Ornstein- 
Uhlenbeck  process  via  eigenfunctions  of  its  associated  Fokker-Plank  equation.  When 
weighted  by  27r\/l  —  the  left  hand  side  of  (3)  becomes  the  probability  density  function 
of  a  bivariate  Gaussian  process,  with  mean  zero,  variances  one,  and  covariance  r.  De¬ 
note  this  density  f{Yi,Y2)-  This  means  that  the  correlation  function  of  the  target  process 
in  the  MNLT  formulation  can  be  written  in  terms  of  that  of  the  input  (Gaussian  pro¬ 
cess)  and  Hermite  polynomials.  Let  rinput  be  the  Gaussian  correlation  function,  so  that 
Tinput{ti-,t2)  =  and  let  the  output  correlation  be  routput{ti-,t2)-  In  particular, 

using  this  and  the  definition  of  correlation  function,  it  can  be  shown  that  for  a  Gaussian 
input  process. 


/OO  poo 

/  h{yi)v{y2)f{Yi,Y2)iyi^y-2)dyidy2 

-OO  j  —OO 
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n 

0 

E{X2) 

1 

V2E{XiX2) 

2 

2  {E{XfX2)  -  E(X2)) 

3 

2V2  (e(a:3x2)  -  3E(a:iA:2)) 

Table  1 :  The  first  few  coefficients  expressed  as  expectations. 


1  Dnpiit  (^1 )  ^2) 

TT  2”n! 

n=0 


n  /  poo 


e  ^  Hn{x)r]{V^x)dx 


=  E 

n=0 

with  the  coefficients  On  given  by 


Dnpiit  (^1 )  ^2)  2 


(4) 


Q'T).  - 


/oo 

e~^  Hn{x)rj{y/2x)dx 

-00 


L 


—00  27r 


e  "2  )  7]{x)dx 


=  E 


(5) 


the  mean  in  (5)  being  with  respect  to  the  Gaussian  marginal  distribution  Xi  =  Ai(0, 1). 
Note  that  if  we  assume  that  the  Gaussian  input  process  is  wide  sense  stationary,  the  output 
process  will  preserve  this  feature,  in  view  of  (4).  Consequently  we  can  write,  in  terms  of 
discrete  time  k  gN, 


OO 

output  {k')  ^  ^ 

n=0 


{j"  inputikif) 


(6) 


Using  the  definition  of  the  Hermite  polynomials,  as  well  as  their  recurrence  relation,  it  can 
be  shown  that  the  first  few  Hermite  polynomials  are  =  1)  Hi (t)  =  2t,  H2{t)  =  —  2 

and  H3{t)  =  8t^  —  12t.  Applying  these  results  to  (5),  we  can  show  the  hrst  few  coefficients 
are  as  given  in  Table  1. 


The  expression  (6)  expands  the  output  correlation  function  in  terms  of  the  Gaussian  input 
correlation  function  as  a  power  series.  It  is  reported  in  [1]  that  this  series  is  rapidly 
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Figure  2:  Relationships  between  input  and  output  autocovariance  coefficients  for  the  Log- 
Normal  distribution. 

convergent  for  a  number  of  target  distributions  of  interest.  This  mapping  between  the  two 
correlation  functions  shows  how,  in  theory,  we  can  control  the  output  correlation  through 
the  Gaussian  input  correlation.  In  essence,  we  require  an  inversion  of  (6). 

As  can  be  seen  in  figure  2,  this  mapping  is  very  much  nonlinear,  and  thus  inversion  in 
the  general  case  is  clearly  not  a  simple  exercise.  Despite  this  difficulty,  as  pointed  out  in 
[1]  and  also  indicated  in  [6],  we  can  establish  some  “rules  of  thumb”  to  determine  which 
Gaussian  process  to  use  as  input.  This  is  addressed  in  the  following  section. 


3.2  Methods  of  Inversion 

Rules  of  thumb  are  now  outlined  for  inversion  of  the  autocorrelation  relationship  derived 
in  the  previous  Section.  Firstly,  it  has  been  reported  in  [1,  6]  that  for  some  clutter  models 
of  interest,  there  is  an  approximate  invariance  between  input  and  output  correlations.  The 
Weibull  distribution  is  a  good  example  of  this.  It  has  been  found  in  many  situations,  there 
is  an  almost  linear  association  between  the  correlation  functions.  This  means  we  can  apply 
a  truncation  of  (6)  to  produce  the  linear  estimate 


routputik)  -  mX2))^  +  {E{XiX2))\nput{k),  (7) 
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and  inversion  is  immediate,  noting  it  may  be  necessary  to  scale  the  linearisation  to  pro¬ 
duce  a  proper  correlation  function  for  the  Gaussian  (which  means  at  /c  =  0  it  must  equal 
I,  the  Gaussian’s  variance).  A  second  approach  is  to  use  ad  hoc  adjustment  of  the  input 
correlation  to  produce  an  output  with  local  dependency  as  described  by  a  sample  corre¬ 
lation.  This  assumes  we  are  not  overly  concerned  about  the  specific  form  of  the  output 
correlation,  but  require  a  correlated  sequence  with  the  correct  marginal  distributions.  A 
third  approach  is  to  use  a  higher  order  polynomial  approximation  to  (6),  such  as  a  cubic  or 
quartic  approximation,  and  apply  inversion  numerically,  scaling  the  appropriate  solution 
to  produce  an  approximate  Gaussian  correlation  function. 


There  is  a  fourth  approach  to  designing  the  appropriate  Gaussian  input  distribution.  The 
Lagrange  Inversion  Theorem  can  be  applied  to  invert  the  implicit  function  defined  through 
the  series  (6).  Suppose  we  have  an  analytic  function  /,  and  that  at  a  point  a  the  derivative 
/'(o)  7^  0.  Then,  the  Theorem  states  that  the  equation  f{w)  =  z  can  be  solved  to  give 
w  =  g{z),  where  g  has  a  series  expansion 


g{z)  =a  + 


OO 

E. 

m=l 


lim  D 


m—1 


w  —  a 


f{w)  -  f{a) 


{z-  f{a)y 


m\ 


(8) 


and  is  the  mth  derivative  with  respect  to  w.  This  can  be  applied  to  (6)  by  selecting 

_  2 

a  =  0  and  letting  f{w)  =  where  bn  =  2^5  noting  that  /  is  analytic  and 

/'(O)  =  ]E2(XiX2)  /  0  in  general.  Here  /(a)  =  /(O)  =  bo.  As  an  example,  four  terms  were 
manually  calculated,  resulting  in  the  approximate  solution 


1  ^2 

f'input{k^  ~  output  {k)  60)  output  {k^  &o) 


(9) 

fh  3  fbi  5b2b3  5bl\  4 

^^4  ^5  J  voutputyk)  bo)  ^  b'^  J  ^0)  • 

Given  a  desired  distribution  Fxj  with  correlation  function  r output,  equation  (9)  gives  an 
indication  of  which  Gaussian  distribution  Xi  should  be  used  to  generate  the  desired  cor¬ 
related  sample.  Observe  that  bo  =  W?{X2)  and  so  the  expansion  (9)  is  in  terms  of  the 
process  autocovariance  function,  assuming  we  have  a  stationary  target  distribution. 


Of  the  four  approaches  outlined  here,  we  will  focus  exclusively  on  the  third,  which  involves 
an  inversion  of  a  partial  series  of  the  correlation  function.  It  is  expected  to  give  more 
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accurate  results  than  a  linearisation,  while  the  approach  based  upon  Lagrange  Inversion 
is  basically  equivalent,  but  requires  more  analytical  work. 

3.3  Lognormal  Distributed  Clutter:  A  Special  Case 

Before  proceeding  to  the  approximate  solutions  to  the  inversion  problem  outlined  in  the 
previous  Section,  we  examine  a  case  where  the  relationship  between  input  and  output 
correlation  can  be  given  in  closed  form.  In  [I]  there  are  several  such  examples  given,  but 
these  are  not  entirely  convincing.  The  example  presented  here  is  perhaps  better  because 
it  indicates  how  the  series  in  (6)  can  be  evaluated  for  a  distribution  of  interest  in  radar. 
We  remark  at  the  outset  that  the  following  simplification  of  (6)  can  be  derived  using 
conditional  probability  without  reference  to  the  series  expansion  (6). 

A  positive  random  variable  X  is  said  to  be  Lognormally  distributed  if  its  logarithm  is  Nor¬ 
mally  distributed.  We  write  X  =  LN{fj,,a^)  where  log(A)  =  This  distribution 

has  been  investigated  in  the  past  as  a  model  for  radar  clutter  [2,  3]. 


Let  Xi  =  Ai(0, 1)  and  X2  =  Then  we  can  select  X2  =  ,  and  con¬ 


sequently  r]{z)  =  Simulating  correlated  returns  is  quite  easy  here  because  of  the 


straightforward  relationship  between  the  input  and  output  distributions.  We  now  evaluate 
the  coefficients  a^;  note  that  by  using  the  form  of  rj  and  applying  a  change  of  variables. 


(10) 


where  W  =  N{a,  1).  We  now  apply  Lemma  3.1,  simplifying  (10)  to 

an  =  a^. 


(11) 


Substituting  into  (6),  we  see  that 


.2 


n=0 
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_  ^(^^rinput{k)+2ii+a^  1-^2) 

where  the  Taylor  series  expansion  for  the  exponential  function  has  been  used.  By  taking 
logarithms  of  (12)  we  can  obtain  the  required  input  correlation.  This  result  is  a  useful  test 
mechanism  for  other  approximate  schemes. 


3.4  Hilbert  Spaces  and  Parseval’s  Identity 

In  the  analysis  to  follow,  we  will  be  exclusively  using  truncation  of  the  series  (6)  and 
applying  an  approximate  inversion  to  it.  It  will  thus  be  important  to  quantify  the  error 
incurred  by  truncation  of  the  series.  In  order  to  achieve  this,  it  is  useful  to  review  some 
Hilbert  space  theory,  which  provides  some  useful  tools  for  this  purpose.  In  particular, 
we  will  require  Parseval’s  Identity  as  well  as  an  understanding  of  the  representation  of  a 
function  in  terms  of  basis  elements. 

Let  X  be  a  vector  space  over  either  the  real  or  complex  numbers,  with  scalar  field  F.  Then 
an  inner  product  [16]  is  a  function 

{■,■):  X  X  X  ^  F 

such  that  all  of  the  following  hold,  for  x,y,z  €  X  and  a  €  F: 


{x  +  y,z) 

=  {x,z)  +  {y,z) 

(13) 

{ax,y) 

=  a{x,y) 

(14) 

{x,y) 

=  {y,x) 

(15) 

{x,x) 

>  0 

(16) 

0  =  {x,  x) 

X  =  0. 

(17) 

That  is,  it  is  a  function  from  a  pair  of  vectors  to  a  scalar  that  is  linear  in  its  first  argument 
(13,14),  conjugate-symmetric  (15),  and  positive-definite  (16,17).  While  (16)  might  look 
somewhat  suspicious  at  hrst,  since  the  inner  product  is  a  complex- valued  function,  the 
application  of  (15)  will  show  that  {x,x)  must  be  real. 
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A  space  X  with  an  associated  inner  product  (•,  •)  is  called  an  inner  product  space.  If  it  is 
also  complete,  it  is  called  a  Hilbert  space.  This  means  every  Cauchy  sequence  in  the  space 
converges  [16],  with  respect  to  the  norm  induced  by  the  inner  product,  to  a  point  in  the 
space.  The  associated  norm  is  ||xp  =  (x,x),  and  is  a  mapping  ||  -  I  :  A  ^  M. 

If  two  vectors  have  an  inner  product  of  zero,  we  call  them  orthogonal.  If  the  elements 
of  a  set  of  vectors  are  orthogonal  to  each  other  element,  we  say  that  the  set  is  mutually 
orthogonal.  If  each  element  has  a  norm  of  one,  we  describe  the  set  as  orthonormal.  If 
this  set  is  sufficiently  large  that  we  can  write  any  vector  as  a  linear  combination  of  its 
elements,  we  call  it  an  orthonormal  basis.  The  representation  in  terms  of  any  particular 
basis  is  unique. 

One  might  then  ask  how  to  find  the  coefficients  of  each  member  of  the  basis  that  are  used 
when  forming  a  particular  vector.  When  the  basis  {e„}  is  orthonormal,  this  is 

X  =  '^en{x,en).  (18) 

n 

Having  found  this  representation,  we  may  then  introduce  the  Parseval  Identity,  which  will 
play  an  important  role  in  determining  the  error  bound  on  our  truncated  correlation  series. 
Since  ||xp  =  {x,x),  it  follows  by  (18)  that 


=  {en{x,en),x) 

n 


=  '^{x,en){en,x) 

n 

=  '^{x,en){x,en) 

n 

=  .  (19) 

n 

In  Euclidean  space,  this  statement  essentially  states  that  even  if  we  select  a  different 
coordinate  system  to  represent  our  vectors,  the  square  of  the  length  of  a  vector  will  still  be 
given  by  the  sum  of  the  squares  of  its  components  as  long  as  the  new  basis  is  orthonormal. 

The  family  of  spaces  that  we  require  for  our  application  is  L‘^{—oo,  +oo),  and  an  appro¬ 
priate  choice  of  basis  functions  is  the  Hermite  polynomials  with  respect  to  the  Gaussian 
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X 


Figure  3:  The  zeroth,  first,  and  sixth  Hermite  functions. 
weight  function  (see  figure  3),  given  by 

Gnix)  =  .  \  e~^'^^Hn{x),  (20) 

where  the  inner  product  on  this  space  is  given  by 

/  +  CXD  _ 

fix)g{x)  dx  (21) 

-OO 

(see  [16]  for  details  of  this,  including  justification  that  the  sequence  {en{x)}  forms  an 
orthonormal  basis). 

We  shall  see  that  the  effect  of  the  MNLT  is  described  most  simply  when  it  is  considered 
in  terms  of  these  functions. 

Before  doing  so,  however,  we  return  once  more  to  Parseval’s  theorem,  and  its  interpretation 
in  the  context  of  functions  in  for  a  general  interval  I.  The  norm  of  a  function  /  is 

ll/p  = 

=  I  l/(t)P  dt  (22) 

The  square  of  the  norm  of  a  function  of  is  therefore  the  total  energy  of  the  function. 

Expanding  (19)  using  (21,  22),  we  see  that  (where  {e^}  is  an  arbitrary  orthonormal  basis) 

[\fitfdt  =  j;i(/.e„>p.  (23) 

n 

We  can  write  the  coefficients  for  each  element  of  the  basis  as  a  discrete  sequence 

Cn  =  {...,(/,  e_i),  (/,  ep),  (/,ei),...}, 
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and  therefore  Parseval’s  theorem  is  the  statement  that  the  energy  in  this  discrete  signal  is 
equal  to  the  energy  in  our  original  function. 


3.5  Inner  Product  Representation 


The  purpose  here  is  to  show  how  the  correlation  series  (6)  can  be  expressed  in  terms  of  the 
inner  product  and  associated  basis  functions  of  the  associated  Hilbert  space  L^(— oo,  +oo). 
This  allows  the  application  of  Parseval’s  Identity  to  the  series,  and  perhaps  more  impor¬ 
tantly  grants  an  understanding  of  the  distribution  of  output  autocorrelation  amongst  the 
infinitude  of  terms. 

Note  that  the  coefficients  of  (5)  can  be  written 


The  Hermite  functions  are  given  by 


en{x) 


"t7,! 


-.e  2 


n'.y/'K 


and  we  define 


ii{x) 


e  ^  rj 


allowing  us  to  then  simplify  (24) 


(25) 


and  thus  write  (6)  as 

OO 

f'outik)  —  'y  ^  input • 
n=0 
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Knowing  that  both  ■0  and  are  real,  this  is  equivalent  to 

OO 

rout{k)  =  y^Jinput{kr\{'ilJ,en)\‘^  ,  (26) 

n=0 

thus  putting  the  series  into  a  form  amenable  to  application  of  Parseval’s  Identity. 

As  the  input  Xi  =  A^(0, 1),  we  have  that  rinput{0)  =  IE  (-^i)  =  1-  Similarly,  by  the 
definition  of  autocorrelation,  rout{0)  =  E  Thus,  we  can  write  a  special  case  of  (26): 

CX) 

E(X|)  =  j;|(0,e„)|^  (27) 

Thus  it  follows  from  Parseval’s  Identity  that 

||0||2  =  E(x|),  (28) 

as  the  Hermite  functions  form  a  complete  orthonormal  basis  for  L^(— oo,+oo).  But,  IV’P 
is  the  energy  of  the  signal  0,  and  therefore  it  has  total  energy  equal  to  the  instantaneous 
power  of  the  target  process. 

One  might  perhaps  choose  to  interpret  the  output  autocorrelation  as  the  energy  in  a 
filtered  version  of  0.  In  essence,  the  multiplication  by  x""  corresponds  to  an  attenuation 
by  a  factor  of  x~'^  of  each  term  n.  When  we  attempt  to  bound  the  remainder,  we  do  so 
by  considering  the  power  in  0  that  the  truncated  series  has  failed  to  include. 


3.6  Bounding  the  Truncation  Error 

The  power  series  described  in  (26)  can  be  truncated  and  used  to  determine  the  autocor¬ 
relation  characteristics  of  the  transformed  process.  We  then  desire  to  know  the  error  that 
such  an  approximation  introduces.  Let  gk{x)  be  the  truncation  of  the  series  to  order  k; 
that  is,  let 

k 

5fc(a;)  =  |(0,en)P  ,  (29) 

n=0 

and  the  corresponding  remainder  function  be 

CXD 

rk{x)  =  ^  x"|(0,e„)|^  (30) 

n=k-\-l 
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It  should  be  noted  that  being  orthonormal  implies  that  setting  'ifj  =  Cp  will  result 

in  =  dnp,  and  thus  we  may  produce  a  relationship  between  input  and  output 

autocorrelation  described  by  an  arbitrary  polynomial  by  selecting  an  appropriate  linear 
combination  of  these  functions.  Adding  to  ■0  the  function  ce^  will  therefore  produce  a 
term  with  coefficient  c  at  n  =  k,  demonstrating  that  one  cannot  bound  the  remainder 
without  some  knowledge  of  the  target  distribution. 


For  an  input  autocorrelation  of  magnitude  at  most  one,  as  in  the  case  of  a  standard  normal 
distribution,  the  worst-case  error  occurs  when  all  the  remaining  energy  of  0  is  in  the  next 
(and  least  attenuated)  term.  Irrespective  of  the  distribution  of  energy  in  0,  for  n  >  k  we 
have  that  \x\^  <  and  therefore  can  bound  the  remainder  series: 


rk{x)\ 


n=k-\-l 


oo 

<  |x"||(0,e„)|2 

n=k-\-l 


< 


E 

n=k-\-l 


l(0,e„ 


2 


oo 

n=k-\-l 


(31) 


By  (27),  we  have  that 

9fc(l)  +  ^fc(l)  =  IE  (2f|)  , 

and  thus  that 

\rk{x)\  <  |x|'=+i(E(X|)  -fffc(l)).  (32) 

When  0  is  a  linear  combination  of  Hermite  functions  of  order  at  most  fc  -|-  1,  (0,  e,i)  =  0 
for  n  >  /c  +  2,  causing  all  but  the  k  +  1*^  term  to  vanish,  and  equality  is  achieved: 

rk{x)  =  |a;|^+Vfc(l) 

=  \x\^+^{¥.{xl)-gk{l)). 
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Note  that  as  gkiX)  increases  monotonically  with  k,  as  its  coefficients  are  all  positive. 
Therefore,  E(X2)  —  5fc(l)  will  be  decreasing,  and  thus  we  can  loosen  (32)  slightly  and 
write  (where  n  <  k) 

rk{x)  <  (E(X|)  -5fc(l)) 


<  \xf+^  (E(X|)-5n(l)) 


We  know,  however,  that  go{k)  =  E(W2)^,  and  therefore  that 


rk{x)  <  |x|'=+1(E(X|)-E(X2)2) 


x|^+^Var(X2), 


(33) 


and  thus,  in  order  to  achieve  a  remainder  of  at  most  rk{x),  we  must  have 


^  ^  log  {rk{x)/\av{X2))  _  ^ 
~  log  X 


(34) 


Hence  we  have  rules  of  thumb  enabling  the  determination  of  errors  associated  with  the 
series  truncation. 


3.7  Numerical  Error 


The  bound  given  in  (31)  is  adequate  only  if  we  are  to  discount  errors  introduced  by 
numerical  integration.  While  in  many  cases  this  might  be  sufficient,  we  would  do  well  to 
be  able  to  predict  the  error  introduced  by  our  implementation. 

Let  the  estimated  value  of  each  coefficient  be  denoted  6*,  with  relative  error  at  worst 
That  is,  suppose  |6*  —  |('0,e„)||  <  £r\bn\-  We  have  therefore  an  error  that  is  at  worst  (by 
the  triangle  inequality) 

k 

iA(x)i  < 

k 

n=0 


6*^ 
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k 

<  J.j;|xr|i.;||2C-K-|W.e„)|)| 

n=0 

k 

n=0 

=  (e^  +  2er)  gk  (Ixj) .  (35) 

Noting  that  coefficient  errors  must  be  applied  to  (32),  we  add  the  truncation  error  to  (35) 
and  find  that  gk{Ti[k\)  differs  from  ro[k\  by  at  most 

|e(x)|  <  [E  (X|)  -  (1  -  Srfgkil)]  +  (e^  +  2e^)  gu  (|x|) .  (36) 

In  the  worst-case  of  x  =  1,  this  reduces  to 

|e(l)|  <E(X|)  -(2e2  +  l)5,(l) 

and  thus  for  a  given  worst-case  error  we  have  the  requirement  that 


Er  < 


1  /E(X|)-|e(l)| 

2  V  9k{l) 


(37) 
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4  Implementation 


The  implementation  of  the  MNLT  is  now  addressed.  It  has  been  implemented  as  a  Matlab 
toolbox  using  a  truncated  power  series  for  autocorrelation  control.  This  task  can  be 
split  into  five  main  segments:  the  determination  of  r]{x),  the  calculation  of  the  input 
autocorrelation  function,  production  of  the  correlated  Gaussian  series,  translation  to  the 
target  distribution,  and  validation  of  the  output.  We  address  each  of  these  individually  in 
the  following  subsections. 

Beyond  this,  the  toolbox  contains  integrated  documentation,  and  examples  of  its  use  are 
given  in  Appendix  A. 


4.1  Constructing  the  Transformation 

The  MNLT  function  r]{x)  =  F^^{Fxi{x)),  defined  in  Lemma  2.4,  transforms  a  standard 
Gaussian  random  variable  into  one  of  the  target  distribution.  However,  it  is  clear  that 
Fx2  is  not  analytically  invertible  in  general.  We  thus  require  a  numerical  approximation 
of  its  inverse.  This  was  done  by  sampling  Fx2  uniformly  over  a  bijective  interval  around 
Fx2  =0.5.  This  then  produced  an  approximation  of  the  inverse  using  the  method  of  cubic 
splines  (see  figure  4).  The  implementation  is  shown  in  Section  B.l. 


F{x) 


X 


X 


Figure  Interpolation  is  used  to  determine  values  of  the  inverse  CDF  between  samples. 


UNCLASSIFIED 


23 


DSTO-TR-2517 


UNCLASSIFIED 


Ttarget 


A 


Per- Realisation 


Initialisation 


Figure  5:  The  means  by  which  the  implementation  correlates  the  output. 

4.2  Determining  the  Input  Autocorrelation 


We  could  use  (5)  to  compute  the  values  of  \{Tp,  e„)p,  using  the  Monte  Carlo  approximation 

.  N  .  s 

l(V’,en)|2  ^  ,  (38) 

n=l  ^  ^ 

where  Xn  are  realisations  of  the  standard  normal  distribution.  This  was  used  at  first, 
however  performance  and  accuracy  were  unsatisfactory. 


Alternatively,  we  may  use  a  deterministic  method  of  numerical  integration.  The  integral 
form  of  (5)  can  be  written  as 

^+oo  ^ 

/  e  ^  f{x)  dx. 

J  —  OO 

Such  an  integral  can  be  estimated  numerically  via 


'^^kfixk)  +  R{f), 


k=l 


where  x^  are  the  roots  of  the  Hermite  polynomial,  and 


Au  — 


(39) 


(40) 


[H'n{xk)f  ■ 

(See  [24]  for  details  of  this).  This  formula  is  exact  for  the  case  when  f{x)  is  a  polynomial 
of  order  at  most  2n  —  1.  The  remainder  of  this  approximation  is  given  by 

n!-/7r 


R{f)  = 


(2n) 
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where  c  is  the  point  at  which  takes  its  maximum. 

A  22-point  numerical  integration  using  this  method  significantly  outperformed  the  two- 
million  point  Monte  Carlo  simulation.  However,  difficulties  in  calculating  (40)  when  using 
a  large  number  of  points  prevented  the  use  of  this  approximation  when  high  accuracies 
were  required. 

Hence  an  adaptive  technique  was  examined;  specihcally,  the  Gauss-Kronrod  algorithm  [22], 
whose  implementation  is  included  with  Matlab.  This  provided  a  high  degree  of  accuracy, 
and  is  thus  used  in  our  implementation. 

Having  done  this,  we  now  have  a  finite  polynomial 

Qkif  input)  =  \{'4’ 1  ^o)\  +  Dnpiitl  (V’)  Cl)  I  +  (V’l  ^2)  j  +  '  '  '  +  Gnpiil  I  Cfc  )  j 

which  can  be  solved  numerically. 

Our  next  task  in  determining  the  autocorrelation  of  the  input  sequence  is  to  find  the  desired 
autocorrelation  of  the  output  sequence.  As  we  are  using  a  truncation  of  (4),  naively 
converting  the  autocovariance  coefficient  to  autocorrelation  using  the  known  mean  and 
variance  of  the  target  distribution  will  produce  a  non-standard  Gaussian  input.  Whenever 
gk{x)  /  Toutputix),  as  in  almost  all  cases,  we  have  that 

9k{l)  <  Toutputil)  =  E  (X|)  =  Var(X2)  +  E  (Xa)'  (41) 

and  hence  solving  gkix)  for  the  right-hand  side  of  (41)  would  produce  rinput{0)  /  1. 

This  issue  can  be  addressed  by  scaling  the  target  autocovariance  by  a  constant.  Multiplying 
P  by  (5fc(l)  —  E(X2)^)/Var(Al2)  we  instead  find  the  equation 

gkirinput)  =  P  (5fc(l)  -  E  (^2)^^  -h  E  {X2f  , 

where  p  denotes  the  target  autocovariance  coefficient.  In  the  case  when  p  =  1,  it  is  clear 
that  Tinput  =  1  is  indeed  a  valid  solution.  As  p  =  1  by  necessity  at  zero-lag,  we  can  then 
see  that  the  resulting  normal  distribution  will  indeed  be  standard,  and  thus  p(^i)  will 
produce  the  correct  output  distribution. 
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4.3  Simulating  the  Correlated  Process 

In  this  brief  section  we  clarify  how  the  simulation  of  a  correlated  Gaussian  process  can  be 
achieved.  The  key  to  this  is  the  well-known  Cholesky  Decomposition  of  a  covariance  matrix 
[25].  Since  a  covariance  matrix  is  positive  definite,  such  a  decomposition  will  always  exists. 
This  factorisation  provides  a  matrix  A  such  that  if  R  is  the  covariance  matrix,  R  =  AA^, 
where  the  latter  matrix  is  the  transpose  of  A. 

The  key  idea  in  simulating  a  desired  correlated  Gaussian  process  is  the  following.  Gonsider 
the  effect  of  a  linear  transformation  on  a  vector  of  independent  Gaussian  variables.  Let 
X  =  [Xi  •  •  •  be  a  column  vector  of  standard  independent  Gausssian  variables,  A  be  an 
n  X  n  matrix,  and  Y  =  [Ti  •  •  •  l^j  =  AX.  Then  it  is  relatively  simple  to  show  that  Y  is  also 
a  Gaussian  process  with  covariance  matrix  AA^  =  R.  Hence,  once  the  Gholesky  factor 
matrix  A  is  identified,  we  can  simulate  the  correlated  process  using  a  matrix  of  independent 
and  identically  distributed  standard  Gaussian  random  variables  and  A.  Specific  to  our 
simulation  problem,  we  can  generate  a  correlated  vector  Z  in  our  target  distribution  with 

Z  =  r,  (AX) , 

where  r/(Y)  is  taken  to  be  the  vector  [rj{Yi)  ■  ■  ■  rj{Yn)]- 


4.4  Testing 

Testing  and  validation  of  the  simulator  built  to  implement  the  MNLT  was  extensive. 
To  test  the  capability  of  the  simulator  to  generate  uncorrelated  processes,  we  used  the 
Kolmogorov-Smirnov  test  [12,  23]  to  verify  that  the  distribution  of  the  output  matched 
the  theoretical  distribution  at  a  confidence  level  of  a  =  0.01,  using  the  Exponential  and 
KK-Distributions  with  10^  elements  in  each  realisation. 

For  correlated  processes  the  Kolmogorov-Smirnov  test  is  not  an  appropriate  measure  of 
the  matching  of  probability  distributions,  because  it  only  applies  to  realisations  that  are 
independent  and  identically  distributed.  Gonsequently  a  50-bin  goodness  of  fit  test  [12] 
was  used  at  a  confidence  level  of  a  =  0.01.  To  provide  an  adequate  sample,  the  test  was 
run  against  the  concatenation  of  1000  realisations,  each  of  length  4000. 
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These  confidence  levels  are  relatively  low  as  the  tests  are  intended  to  be  run  repeatedly 
with  various  distributions  as  part  of  an  automated  test  suite,  and  the  use  of  a  =  0.05  is 
therefore  likely  to  produce  an  unacceptable  rate  of  false  positives. 
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5  Simulation  Results 


We  are  now  in  a  position  to  examine  the  results  of  the  application  of  the  MNLT  to  a 
number  of  clutter  models  of  interest. 

Recall  that  in  Section  3.3  we  found  a  closed-form  relationship  between  the  input  and 
output  autocorrelations  of  the  MNLT  when  used  to  simulate  the  Log-Normal  distribu¬ 
tion.  Figure  6  provides  a  comparison  of  the  result  of  the  MNLT  with  a  series  truncation 
performed,  and  the  corresponding  analytic  results.  The  first  plot  in  figure  6  shows  the 


t 


Figure  6:  Simulation  of  a  Correlated  Lognormal  Process 

amplitude-squared  simulation  results,  the  second  plot  shows  the  error  in  cumulative  dis¬ 
tribution  function,  the  third  plot  shows  the  autocovariance  function  and  the  fourth  plot 
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shows  the  autocovariance  error.  The  MNLT  uses  a  fourth-order  power  series  truncation. 
The  distribution  used  here  is  the  standard  Lognormal  LN (0,1).  We  see  that  there  is  a 
peak  error  around  zero  in  the  approximate  inversion,  corresponding  to  larger  target  and 
input  autocovariances.  This  error  reached  a  maximum  value  of  approximately  5  x  10“^.  It 
is  clear  here  that  the  MNLT  performs  very  well.  Note  that  we  are  using  an  autocovariance 
coefficient  given  by  the  expression  r{k)  = 


Figure  7:  Simulation  of  a  Correlated  Weibull  Process 

Next  we  examine  the  case  of  generating  correlated  Weibull  clutter;  figure  7  shows  the 
four  relevant  plots.  The  first  shows  the  intensity  measurements,  the  second  shows  the 
cumulative  distribution  function  error,  the  third  shows  the  autocovariance  while  the  fourth 
is  the  corresponding  error.  In  this  simulation,  the  target  distribution  is  the  Weibull  with 
cumulative  distribution  function  F{t)  =  1  —  The  correlation  function  is  r{k)  = 
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.  The  output  in  plots  three  and  four  is  the  result  of  the  MNLT  processing  of  a  fourth- 
order  correlation  function  as  before.  We  observe  that  the  MNLT  is  working  very  well  here 
as  for  the  Lognormal  example  considered  previously. 

The  next  example  we  consider  is  for  simulation  of  correlated  KK-clutter  returns.  Recall 
that  the  KK-Distribution  is  the  weighted  sum  of  two  K-Distributions;  we  can  write  its 
cumulative  distribution  function  as 


FKK{t;k,U,Ci,C2)  =  (1  -  k)FK{t]V,Ci)  +  kFK{t-,V,C2), 


where  0  <  fc  <  1  is  the  distribution  mixing  parameter,  ci  and  C2  are  the  two  scale  pa¬ 
rameters  and  u  is  the  KK-Distribution’s  shape  parameter.  The  individual  K-Distributions 
used  in  the  mixture  have  cumulative  distribution  function  given  by 


FK{t-,iy,Ci)  =  1  - 


{citYKy{cit) 


for  i  e  {1,  2}. 


We  consider  the  case  where  the  KK-Distribution  is  characterised  hy  k  =  0.01,  v  =  3.4119, 
Cl  =  110  and  C2  =  85.  The  correlation  function  is  r{k)  =  cos(0.l7rfc)e“^'^^.  The  first 
simulation  is  in  figure  8,  which  examines  the  results  of  a  linear  approximation  to  the 
autocorrelation  relationship.  The  autocovariance  errors  are  quite  large;  this  is  perhaps  to 
be  expected,  as  it  is  the  higher-order  terms  of  (4)  that  encode  the  distortion  produced  by 
the  MNLT. 


The  use  of  a  quadratic  approximation  as  shown  in  figure  9  improves  dramatically  the 
quality  of  the  simulation,  reducing  the  peak  autocorrelation  error  to  nearly  one  quarter 
that  achieved  through  a  linear  approximation. 

A  fourth-order  approximation  is  shown  in  figure  10,  which  further  decreases  error  to  ap¬ 
proximately  one  half  that  achieved  using  a  quadratic.  It  is  clear  that  increasing  the  quality 
of  the  approximation  will  reduce  the  autocovariance  error.  Figure  11  demonstrates  the 
simulation  of  another  KK-distribution,  showing  that  the  simulator  performs  well  for  a 
wide  range  of  parameters.  In  this  case,  we  chose  k  =  0.3,  v  =  2.5,  ci  =  1  and  C2  =  3,  as 
well  as  an  autocovariance  coefficient  of  r{k)  =  (1  -|-  t)e~^.  We  observe  the  MNLT  again 
works  very  well  for  the  KK-Distribution. 
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Figure  8:  Simulation  of  a  KK  Process  With  Linear  Autocorrelation  Correction 


Calculated  values  of  the  series  coefficients  and  remainders  are  shown  in  Appendix  C;  we 
see  that  for  the  KK-Distribution,  the  series  indeed  converges  quite  rapidly. 
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Figure  9:  Simulation  of  a  KK  Process  With  Quadratic  Autocorrelation  Correction 


UNCLASSIFIED 


33 


DSTO-TR-2517 


UNCLASSIFIED 


0.01 


0  200  400  600  800  1000 


t 


Figure  10:  Simulation  of  a  KK  Process  with  Quartic  Autocorrelation  Control 
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Figure  1 1 :  Simulation  of  another  KK  Process  with  Quartic  Autocorrelation  Control 
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6  Conclusions 

The  Memoryless  Nonlinear  Transform,  which  enables  the  simulation  of  correlated  sea 
clutter  from  a  prescribed  distribution  with  given  autocovariance  function,  has  been  in¬ 
troduced.  This  process  involves  converting  a  correlated  Gaussian  stochastic  process  (the 
input  process)  to  a  desired  process  with  given  autocovariance  function  (the  output  pro¬ 
cess).  A  new  technique,  which  truncates  the  mapping  between  input  and  output  process 
autocovariances  and  then  performs  an  inversion,  has  been  introduced.  This  enables  the 
determination  of  which  input  Gaussian  process  should  be  used  to  generate  the  desired 
output  stochastic  process. 

A  simulator  has  been  constructed  in  Matlab,  and  testing  of  some  examples  showed  excellent 
results.  In  particular,  we  showed  the  MNLT  works  very  well  for  earlier  models  of  sea 
clutter,  namely  Lognormal  and  Weibull  Distribution  models.  It  was  also  demonstrated 
that  the  MNLT  worked  very  well  for  the  generation  of  correlated  KK-Distributed  clutter. 
Other  tests,  including  an  analysis  for  the  case  of  the  K-Distribution,  further  supported 
the  validity  of  this  simulator. 

The  simulator  provides  DSTO  with  a  capability  that  can  be  integrated  into  radar  models 
under  development  in  the  Microwave  Radar  Branch.  This  should  prove  to  be  useful  for 
the  longer  term  goals  of  testing  high  grazing  angle  detection  schemes  operating  under 
correlated  sea  clutter  returns. 
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Appendix  A  Examples  of  Toolbox  Usage 

This  section  describes  some  basic  use  of  the  toolbox.  More  detailed  documentation,  in¬ 
cluding  a  tutorial  and  a  short  explanation  of  the  theory  of  operation,  is  available  in  the 
Matlab  help  browser  and  in  the  help  for  individual  functions  (see  Section  B). 


A.l  Generating  an  Uncorrelated  Series 

This  example  shows  how  one  might  simulate  a  uncorrelated  exponential  process.  Using  the 
expcdf  function  from  the  Matlab  Statistics  Toolbox,  it  produces  one  million  uncorrelated 
samples  distributed  according  to  Exp{l/2).  The  code  in  Listing  1  is  also  available  in  the 
toolbox  as  example/uncorrelated,  m 


Listing  1:  Uncorrelated  Simulation 


%  Generate  an  uncorrelated  exponential 

%  some  features  of  the  process  as  well 

distribution  ,  and  show 

as  a  histogram  . 

ctx  =  Simulator  (@(  t )  expcdf(t,2),  le6  , 

X=  simulate  (  ctx  ) ; 

0.01); 

fprintf(  ’  Output  ^  s  ize  :  \  n  ’  ) ; 

disp ( size  (X) ) ; 

fprintf(  ’  Output  ^mean :  \  n  ’  ) ; 

disp  (mean(X) ) ; 

fprintf(  ’  Output  ^variance  :  \  n  ’  ) ; 

disp  (  var  (X) ) ; 

hist (X,  1000); 
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Listing  2:  Uncorrelated  Simulation  Output 


Output  size  : 

1 

1000000 

Output  mean : 

2.0011 

Output  variance 

1  • 

4.0025 

0  5  10  15  20  25  30 


Figure  Al:  Uncorrelated  Simulation  Output  Histogram 
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A. 2  Generating  a  Correlated  Series 


This  example  shows  how  one  might  simulate  an  correlated  exponential  process.  The 
interface  is  the  same  as  the  uncorrelated  case,  however  memory  and  computational  con¬ 
siderations  restrict  us  to  a  shorter  sequence.  The  code  in  Listing  3  is  also  available  in  the 
toolbox  as  example/ correlated,  m. 


Figure  A2:  Correlated  Simulation  Output 
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Listing  3:  Correlated  Simulation 

%  Generate  a  correlated  exponential  distribution  ,  and  show 
%  some  features  of  the  process  as  well  as  a  histogram  and 
%  plot  . 

r  =@(t)  cos  ( t  *  pi  / 1 0 ) .  *  exp(— t  ) ; 

ctx  =  Simulator  (@(  t )  expcdf  ( t  ,  2 )  ,  r(0:999),  0.01); 

X=  simulate  (  ctx  ) ; 

subplot  (2,1,1); 
hist (X) ; 

t  it  le  (  ’ Output  „ Histogram  ’  ) ; 

%  Average  the  autocovariance  over  Ik  realisations  . 

lags  =  —50:50; 

cov  =  zeros ( size ( lags  )) ; 

for  i  =  1:10  0  0; 

X=  simulate  (  ctx  ) ; 

cov  =  cov  +  xcov(X,  50,  ’coeff’); 

end ; 

cov  =  cov  /  1000; 

subplot  (2  ,1  ,  2  ) ; 

scatter ( lags  ,  cov,  ’.b’); 

hold  on;  plot ( lags  , r (abs ( lags )),’ r ’) ;  hold  off; 
box  on ; 

legend  (  ’  Output  ’  ,  ’Target  ’  ) ; 

xlabel ( ’ Lag ’ ) ; 

ylabel  (  ’  Autocovariance  ’  ) ; 

title  (  ’ Output ^and^Target ^ Autocovariances  ’  ) ; 


autocovariance 
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Appendix  B  Simulator  Toolbox  Code  Listings 

function  obj  =  Simulator  (F ,  r,  varargin) 

%  Simulator  Simulate  an  autocorrelated  time  series. 

% 

%  ctx  =  Simulator(F,  N); 

%  ctx  =  Simulator  (F,  autocovariance  ( 1  :N) )  ; 

% 

%  ctx  =  Simulator(F,  r,  w); 

% 

%  ctx  =  Simulatorf  ,  ’ approximationOrder  n); 

%  ctx  =  Simulator(  ...  ,  ’  approximationTol  n); 

%  ctx  =  Simulator(  ...  ,  ’  controlAutocorrelation  c); 

% 

%  Example  : 

% 

%  %  Generate  uncorrelated  exponential  variables. 


% 

E 

=  @(t)  (  1—  exp(—  t )  )  .*  (t  >=  0); 

% 

ctx 

=  Simulator  (F,  1000); 

% 

% 

XA 

=  simulate  (  ctx  )  ; 

% 

X.2 

=  simulate  (  ctx  ) ; 

% 

% 

%  Generate  correlated  exponential  variables 

% 

E 

=  @(t)  (  l—exp(—t)  )  .*  (t  >=  0); 

% 

r 

=  @(t)  cos  (t*pi  / 10)  .*  exp(—t  / 10) ; 

% 

ctx 

=  Simulator  (F,  r(OAOOO)); 

% 

X 

=  simulate  (  ctx  ) ; 

% 

%A  simulator  object  can  be  used  to  generate  autocorrelated 
%  sequences  of  arbitrary  distributions.  Each  instance 
%  is  immutably  tied  to  a  distribution  ,  sequence  length  and 
%  autocovariance  . 

% 

%  The  simulate  method  produces  a  realisation  of 
%  an  autocorrelated  time  series  of  the  distribution  with  the 
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%  given  CDF. 

% 

%  When  called  with  a  vector  r  of  length  N,  correlated  vectors 
%  of  length  N  will  be  generated  with  the  provided 

%  autocorrelation  coefficients  (ie.  r(l)  should  be  equal  to  one). 

% 

%  When  called  with  a  scalar  N  in  place  of  an  autocovariance  vector, 

%  uncorrelated  vectors  of  length  N  will  be  generated  . 

% 

%  The  w  argument  determines  the  spacing  between  samples  of  the 
%  CDF.  If  unspecified  ,  this  defaults  to  0.01.  Set  this  to  a  smaller 
%  value  for  CDFs  that  change  quickly . 

% 

%  The  approximationOrder  argument  determines  the  order  of 
%  the  polynomial  approximation  used  when  controlling  the 
%  autocovariance  .  If  not  specified,  this  defaults  to  four. 

% 

%  The  approximationTol  argument  specifies  the  relative  tolerance  of 
%  the  autocorrelation  approximation  coefficients. 

%  This  defaults  to  le  —  9. 

% 

%  The  controlAutocovariance  argument  disables  output 
%  autocovariance  control .  If  this  is  false,  then  the  input 
%  autocovariance  will  be  used  directly  with  the  MNLT.  This 
%  might  be  useful  in  cases  where  the  effect  of  the  MNLT  on 
%  the  output  autocovariance  can  be  determined  analytically. 

% 

%  See  also  Simulator .  simulate 

%  Initial  values  to  make  sure  the  fields  are  created . 

obj  .U  =  []  ; 

obj . correlated  =  0; 

obj  .  variance  =  1 ; 

obj  .  approximationOrder  =  0; 
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obj  .  approximationTol  =  0; 
obj  .  controlled  =  false  ; 
obj  .  w  =  0  ; 
obj  .F  =  0; 
obj  .  Finv  =  0  ; 

%  Produce  a  dummy  class  if  called  with  no  arguments  so  that  we 
%  can  load  Simulator  objects  from  a  file  . 

if  nargin  =  0 

obj  =  class(obj,  ’Simulator’); 

return ; 

end 

p  =  inputParserCompat  ( ) ; 

p  =  addRequired  (p  ,  ’F’,  @(x)  isa(x,  ’  function.handle  ’ )  ); 

p  =  addRequired  (p ,  ’r’,  @isnumeric  ); 

p=addOptional(p,  ’w’,  0.01,  @(x)  isnumeric(x)  &&  isscalar  (x)  &fe  x  >  0) 
p  =  addParamValue  (p  ,  ’  control  Autocovariance  ’  ,  1,  @(t)  t==0||t==l); 

p  =  addParamValue  (p ,  ’  approximationOrder  ’  ,  4,  @(t)  0  =rem(t,l)  ); 
p  =  addParamValue  (p ,  ’approximationTol’,  le— 9,  @(t)  0  <  t  &&:  t  <  1); 
p  =  parse (p ,F , r , varargin {:}) ; 

results  =  Results(p); 

obj  .  approximationOrder 
obj  .  approximationTol 
obj  .  controlled 
obj  .w 
obj  .F  =  F; 

%  We  begin  by  inverting  the  CDF. 

[a  b]  =  cdf. bounds  (F ,  obj.w); 
obj  .  Finv  =  inversef(F,  a,  b,  obj.w); 

r.size  =  size (r  ) ; 


=  results  .  approximationOrder  ; 

=  results  .  approximationTol  ; 

=  results  .  controlAutocovariance 
=  results  .w; 
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if  isscalar(r) 

%  We  have  been  given  a  scalar,  and  so  generate 
%  uneorrelated  values.  Since  we  need  no  correlation 
%  matrix,  we  save  memory  and  only  store  the  number  of 
%  values  to  generate  . 

obj  .U  =  r  ; 

obj . correlated  =  0; 

else 

%  We  have  an  auto  correlation  function  /matrix ,  and  so 
%  generate  correlated  data. 

%  If  the  output  is  to  be  correlated  ,  then  we  must 
%  provide  an  determine  the  input  covariance  that 
%  is  needed  to  produce  the  desired  output  correlation  . 
if  obj  .  controlled 

r_input  =  invert-autocorrelation  (r,  obj.Finv,  ... 

obj  .  approximationOrder  , 
obj  .  approximationTol  ) ; 

else 

%  We  do  not  need  to  control  the  autocovariance  . 
r_input  =  r  ; 

end 

obj  .  variance  =  r_input(l); 

obj .correlated  =  1; 
if  isvector(r) 

%  We  have  been  given  a  vector  ,  implying  that  the 
%  process  is  stationary  . 

%  If  the  process  is  stationary  then  we  must  construct 
%  a  correlation  matrix. 

% 
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%  DOC.GAUSSIAN.GENERATE 
C  =  ones  ( length  (  r  )) ; 
for  i  =  1 :( length  (  r  ) ) 

for  j  =  1 :( length  (  r  ) ) 

C(i,j)  =  r.input  (abs(i-j  )  +  1); 

end 

end 

%  We  use  a  Gholesky  factorisation  so  that  by 
%  multiplying  a  vector  of  normal  variables  by  U  we 
%  produce  a  vector  of  Gaussian  variables  with  the  given 
%  correlations  . 

[obj.U  p]  =  chol(C); 

%  Gheck  whether  G  was  positive —definite  . 
if  p  =  1 

error(  ’Simulator:Construct:BadCovariance’ 

[  ’  Invalid^autocovariance  ;  ^  ’  .  .  . 

’  Gaussian^  CO  variance  ^matrix  „must„  ’  ... 

’be^positive^definite.  ’]); 

end 


elseif  r_size(l)  =  r_size(2) 

%  TODO:  Simulate  processes  with  nonstationary  autocovariances . 
error  (  ’  Simulator  :  Construct  :  Notimplemented  ’  ,  ... 

’Non— st  at  ionary  „correlations^not  „  implemented  .  ’  ) ; 

else 

error  (  ’  Simulator  :  Construct  :  NonSquareMatrix  ’  ,  ... 

’Autocorrelation  „must  „be  „scalar„or„vector.  ’); 

end 


end 


obj  =  class(obj,  ’Simulator’); 
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B.l  CDF  Inversion 

function  [a  b]  =  cdf_bounds  (F ,  w) 

%  CDFJ30UNDS  Determine  CDF  bounds  of  inversion . 

% 

%  The  hounds  are  discovered  by  placing  both 
%  at  a  point  near  F~^{0.5) ,  then  decreasing  and 
%  increasing  the  lower  and  upper  bounds  respectively 
%  until  we  reach  a  point  at  which  moving  further  results 
%  in  no  change . 

%  Arbitrary  initial  starting  value. 
a  =  0 ; 

%  We  first  locate  ourselves  at  approximately  the  median. 
while  (F(a+w)  <  0.5) 

n  =  0 ; 

while  F(a+w*2"n)  <  0.5 
n  =  n  +  1 ; 

end 

a  =  a  +  w*2 "  n ; 

end 

while  (F(a)  >  0.5) 
n  =  0 ; 

while  F(a— w*2~n)  >  0.5 
n  =  n  +  1 ; 

end 

a  =  a  —  w*2  ~  n ; 

end 

b  =  a ; 

%  Then,  we  expand  the  upper  and  lower  bounds  outwards  until 
%  either  we  reach  an  invalid  value  or  the  limit  of  numerical 
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%  resolution . 

% 

%  The  goal  here  is  that  in  the  range  over  which  we  invert  F, 

%  the  inverse  is  a  single— valued . 
while  (F(a— w)  >=  0  &&  F(a)  ~=  F(a+w)) 
a  =  a  —  w; 

end 

while  (F(b+w)  <=  1  &&:  F(b)  ~=  F(b+w)) 
b  =  b  +  w; 

end 

%  Some  safety  margin. 
a  =  a  +  w; 
b  =  b  —  w; 

end 

function  f  =  inversef(F,  a,  b,  width) 

%  INVERSEF  Approximate  an  inverse  CDF. 

% 

%  The  inversef  function  approximates  an  inverse  by  sampling 
%  the  function  between  the  bounds  given  at  points  separated 
%  by  the  width  argument.  It  then  interpolates  between 
%  these  points  to  return  a  function  handle  based  on  the 
%  piecewise  polynomial  generated  by  the  interpolator. 

% 

%  The  function  returned  will  he  equal  to  a  at  points  less  than 
%  F(a),  and  b  at  points  greater  than  F(bO. 

% 

%  See  also  interpl. 

X  =  a  :  width  :  b  ; 
y  =  F(x); 

%  Generate  a  piecewise  polynomial  from  our  samples . 
pp  =  interpl (y , X spline  ’pp  ’) ; 
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%  Construct  our  function  handle  . 

f  =  @(t)  ppval(pp,  t).*(  (t  >=  F(a)).*(t  <=  F(b))  )  + 

(t  <  F(a))*a  +  (t  >  F(b))*b; 

end 
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B.2  Correlation  Control 

function  b_n  =  an_coeffs  ( eta  ,  max.n  ,  reltol) 

%  AN.COEFFS  Determine  eoefficients  for  the  autocorrelation  series. 

% 

%  The  an_coeffs  takes  as  input  a  function  handle  eta,  and  returns 
%  the  first  max-n  coefficients  using  an  adaptive  quadrature. 

n  =  f  li  p  Ir  (  0  :  max_n  )  ; 

b_n  =  zeros  ( max_n  +  l ,  I); 
first  =  true  ; 
smallest_place  =  I; 
for  k  =  n; 

H  =  hermite_poly  (k  ) ; 

integrand  =  @(t)  z e r o _i n f i n i t i e s (  ... 

poly val  (H/sqrt ( f ac t or i al ( k) *2 ~ k ) ,  t)  ... 

.*  exp(— 1.~2)  .*  et a ( t * sqrt ( 2  ) ) ) ; 

b_n  (niax_n— k+1)  =  quadgk  ( integrand  ,  — Inf  ,  Inf  ,  ... 

’RelTol’  ,  reltol  )''2/pi; 

end 

function  r_input  =  invert-autocorrelation  ( r  ,  Finv  ,  approximationOrder  ,  ... 

tol ) 

%  INVFRTjiUTOCORRFLATION  Invert  the  I—>0  autocorrelation  relationship  . 

% 

%  The  invert-autocorrelation  function  determines  the  input 
%  autocorrelation  of  the  standard  Gaussian  input  to  produce  the  desired 
%  output  autocorrelation  using  the  method  derived  by 
%  Tough  and  Ward  (1999). 

% 

%  Numerical  integration  is  used  to  determine  the  first  few  co  efficients 
%  of  the  power  series  relating  the  autocorrelations  .  The  provided 
%  autocovariance  co  efficients  are  converted  to  our  desired  output 
%  autocorrelation  ,  and  the  polynomial  made  by  truncating  the  power 
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%  series  is  solved,  yielding  our  output  coefficients. 

autocorrelation_poly  =  an_coeffs  (@(  t )  Finv  (0 . 5 *  ( 1  +  erf  ( t  / sqrt  ( 2  ) ) ) )  , 

approximationOrder  ,  ... 

tol  ) ; 

%  We  have  been  given  autocovariance  coefficients  ,  but  in  fact  need 
%  an  autocorrelation  .  Fortunately  ,  we  can  convert  between 
%  them. 

% 

%  While  E[X_a  X_(a—k)]  =  r[k]  *  var(X)  +  mean(X)''2,  we  must  have 
%  r^input  [0]  =  1  (which  our  approximation  does  not  preserve),  and  thus 
%  scale  r[k]  by  a  constant  to  yield  this  autocorrelation  ,  whose 
%  inversion  at  zero  *is*  equal  to  one. 

mean  .squared  =  autocorrelation_poly(  length  (  autocorrelation.poly  ) ) ; 
r  =  r  *  ( poly  val  (  aut  ocorr  elat  ion.poly  ,1)  ~  mean.squared  )  ... 

+  mean.squared  ; 

%  All  this  done,  we  can  perform  the  inversion  itself. 
r.input  =  zeros  (  size  (  r  )) ; 
for  i  =  l:length(r) 

r_input(i)  =  solvepoly(autocorrelation_poly,  r(i)); 

end 

end 

function  x  =  solvepoly(p,  r) 

%  SOLVEPOLY  Solve  a  polynomial  equated  to  a  constant . 

% 

%  The  solvepoly  function  solves  equations  of  the  form 

% 

%  r  =  a-n  x"n  +  ...  +  a-1  x  +  a-0 

% 

%  returning  the  greatest  real  solution. 

%  This  is  equivalent  to  an  equation  of  the  form 

% 
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%  0  =  a-Ti  x'n  +  ...  +  a_0  —  r 

% 

%  which  is  already  solvable. 
new_p  =  p ; 

new.p  ( length  (p ) )  =  new_p  ( length  ( p ) )  —  r; 

candidate -Solutions  =  roots  ( new_p  ) ; 

%  We  restrict  ourselves  to  real  solutions. 
real-Solutions  =  imag(  candidate  .solutions  )  =  0; 
candidate  .solutions  =  candidate-solutions  (  real-solutions  ) 

%  We  select  the  greatest  of  these  solutions. 

X  =  max(  candidate  .solutions  ) ; 

end 


56 


UNCLASSIFIED 


UNCLASSIFIED 


DSTO-TR-2517 


B.3  Series  Generation 


function  CorrelatedN  =  generate  _gaussian  (  obj  ) 

%  GENERATE^GAUSSIAN  Generate  correlated  Gaussian  variables  . 

% 

%  The  g enerate -g aus sian  method  produces  a  vector  of  correlated 
%  standard  Gaussian  variables  based  upon  the  correlation 
%  characteristics  specified  during  initialisation. 

% 

%  The  vector  generated  will  be  of  the  size  of  the  correlation 
%  sequence  ,  the  rank  of  the  correlation  matrix,  or  the  size 
%  given  during  initialisation  ,  as  appropriate  . 

% 

%  Example  : 

% 

%  ctx  =  Simulator .OldStyle  (E ,  r,  w); 

%  X  =  generatc-gaussian  ( ctx ) ; 

if  obj . correlated  =  0 

CorrelatedN  =  randn(l,  obj.U); 

else 

normal-vars  =  randn(l,  length  (  obj  .U) )  ; 

CorrelatedN  =  normaUvars  *  obj.U; 

end 

end 


function  [CorrelatedU  CorrelatedN]  =  generate_uniform  ( obj  ) 

%  GENERATEJJNIFORM  Generate  correlated  uniform  variables. 

% 

%  The  generate-uniform  method  produces  a  vector  of  correlated 
%  uniform  variables  on  [0,1]  by  manipulating  a  vector  of 
%  correlated  normal  variables. 

% 

%  Example  : 

% 

%  ctx  =  Simulator .OldStyle  (E ,  r,  w); 
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%  X  =  generate-uniform(ctx); 

CorrelatedN  =  generate.gaussian  (  obj  ) ; 

%  CorrelatedU  =  normcdf  (  CorrelatedN ) ; 

CorrelatedU  =  0 . 5  *  erfc  (  CorrelatedN/ sqrt  ( 2  )) ; 

end 

function  [X  G  U]  =  simulate  (  obj  ) 

%  SIMULATE  Generate  correlated  realisation  . 

% 

%  The  simulate  method  produces  a  realisation  of 
%  an  autocorrelated  vector  of  the  distribution  with  the 
%  given  CDF. 

% 

%  Before  calling  simulate  ,  the  CDF  must  be  set  with  the 
%  ’F’  attribute  . 

% 

%  Example  : 

% 

%  F  =  @(t)  (  l—exp(—t)  )  .*  (t  >=  0); 

%  ctx  =  Simulator .OldStyle  (F ,  1000,  0.01); 

%  X  =  simulate  (  ctx  ) ; 

% 

%  In  addition  to  the  realisation  of  the  desired  distribution  , 

%  the  method  returns  the  corresponding  normal  and  uniform 
%  vectors  that  were  used  to  generate  the  distribution. 

% 

%  [X  gaussian  uniform]  =  ctx  .  simulate  () ; 

% 

%  The  relationship  between  these  is  such  that 

% 

%  eta ( gaussian )  =X 

% 

%  and  thus  this  is  sufficient  to  tune  the  input  autocovariance  . 
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%  We  proceed  by  applying  the  inverse  CDF  to  the  uniform 
%  distribution  . 

[U  G]  =  generate_uniform  (  obj  ) ; 

%  We  must  take  the  real  part,  as  otherwise  we  sometimes  get 
%  X  +  Oj  ,  which  Octave  does  not  like. 

X  =  real  (  obj  .  Finv  (U) )  ; 

end 

B.4  Miscellaneous 

function  [H  Hnml]  =  hermite  _poly  (N) 

%  HERMITEJ’OLY  Return  the  n’th  Hermite  Polynomial. 

% 

%  The  hermitc-poly  function  returns  the  coefficients  of  the  n’th 
%  Hermite  Polynomial  as  the  first  value,  and  the  coefficients 
%  of  the  n—l’th  Hermite  Polynomial,  with  a  leading  zero. 

% 

%  The  nonzero  coefficient  of  greatest  order  for  the  polynomial  of 
%  order  n  is  2~n. 

if  N  =  0 

H  =  [1]; 

Hnml  =  [  ]  ; 
elseif  N  =  1 

H  =  [2  0]; 

Hnml  =  [  1  ]  ; 

else 

Hnml  =  [  1  ]  ; 

H  =  [2  0]; 

for  i  =  1:(N— 1) 
old.H  =  H; 

H=  [H*2  0]  -  [0  n  Hnml  *2* i  ]  ; 

Hnml  =  old.H  ; 

end 
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Hnml  =  [  0  Hnml  ]  ; 

end 

end 


B.5  InputParser  Replacement 

The  Matlab  inputParser  class  was  required,  but  this  was  not  available  in  Octave.  We  have 
written  our  own,  included  here. 

function  obj  =  inputParserCompat  ( ip  ) 
if  nargin  =  0 

obj.  required  =  {}; 
obj  .  optional  =  {  }  ; 
obj  .  params  =  {  }  ; 
obj  .  Results  =  {  } ; 

obj  =  class  (obj,  ’inputParserCompat’); 
elseif  strcmp  (  class  ( ip  )  ,  ’inputParserCompat’) 
obj  =  ip  ; 

else 

error  (  ’  inputParserCompat  :  ^ Expect ing^no^ arguments  .  ’  ) ; 

end 

end 

function  obj  =  addRequired  ( this  ,  name,  validator) 

obj  =  this  ; 

param  .  name  =  name ; 

param  .  validator  =  validator; 

id  =  length  (  obj  .  required  )  +  l; 
obj  .  required  {  id  }  =  param; 


end 

function  obj  =  addOptional  ( this  ,  name,  default  ,  validator) 
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obj  =  this  ; 

param  .  name  =  name ; 
param  .  default  =  default; 
param  .  validator  =  validator; 

id  =  length  (  obj  .  optional  )  +  I; 
obj  .  optional  {  id  }  =  param; 


end 

function  obj  =  addParamValue  ( this  ,  name,  default  ,  validator) 

obj  =  this  ; 

param  .  name  =  name ; 
param  .  default  =  default; 
param  .  validator  =  validator; 

id  =  length  (  obj  .  params)  +  l; 
obj  .  paramsj  id  }  =  param; 


end 

function  obj  =  parse(this  ,  varargin) 
obj  =  this  ; 

if  nargin  <  length  (  obj  .  required )  +  l 
error  (  ’ Not  „enough„ parameters  .  ’  ) ; 

end 

%  We  start  by  dealing  with  all  of  the  required  arguments . 
i  =  1; 

for  param  =  obj  .  required 
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value  =  vararginji}; 
param  =  param{l}; 

if  is  a  (param  .  validator  ,  ’function_handle  ’) 
if  'param  .  validator  (  value  ) 

error  ( [  ’  Parameter^  ’  ’  ’  param  .  name  ’  ’  ’..failed^validation  .  ’]) 

end 

end 

obj  .  Results  =  setfield  (  obj  .  Results  ,  param. name,  value); 
i  =  i  +  1 ; 

end 

%  Now  set  the  optional  arguments  that  were  supplied . 
parameter  =  1; 

while  i  <=  length  (  varargin )  &fe  parameter  <=  length  (  obj  .  optional ) 
value  =  vararginji}; 
param  =  obj  .  optional  {parameter  }  ; 

if  is  a  (param  .  validator  ,  ’function_handle  ’) 
if  'param  .  validator  (  value  ) 

error  ( [  ’  Parameter  „  ’  ’  ’  param  .  name  ’  ’  ’^failed^validation  .  ’]) 

end 

end 

obj  .  Results  =  setfield  (  obj  .  Results  ,  param. name,  value); 
i  =  i  +  1 ; 

parameter  =  parameter  +  1; 

end 

%  Set  the  other  optional  arguments  to  their  default  values. 
while  parameter  <=  length  (  obj  .  optional ) 
param  =  obj  .  optional  {parameter  }  ; 

obj  .  Results  =  setfield  (  obj  .  Results  ,  param.  name,  param .  default  ) ; 
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parameter  =  parameter  +  1; 

end 

%  Set  all  of  the  named  arguments  to  their  default  values. 
for  param  =  obj  .  params 
param  =  paramjl}; 

obj  .  Results  =  setfield  (  obj  .  Results  ,  param.  name,  param .  default  ) ; 

end 

%  Now  look  for  the  remaining  named  arguments  . 
while  i  <  length  (  varargin ) 
paramName  =  vararginji}; 
if  ~  ischar  (paramName) 

error  (  ’  Parameter^name^must^beu^au^string  .  ’  ) ; 

end 

matched  =  0; 
for  param  =  obj  .  params 
param  =  paramjl}; 
if  strcmp  (paramName,  param.  name) 

value  =  vararginj  i +1}; 

if  isa  (param  .  validator  ,  ’  function_handle  ’ ) 

if  'param  .  validator  (  value ) 

error  ([’ Parameter^ ’’  ’  param  .  name  ... 

’  ’  ’^failed^validation.  ’]); 

end 

end 

obj  .  Results  =  setfield  (  obj  .  Results  ,  param. name,  value); 
matched  =  1; 


UNCLASSIFIED 


63 


DSTO-TR-2517 


UNCLASSIFIED 


break ; 

end 

end 

if  matched  =  0 

error  ([’’’  ’  paramName  ’  ’  ’  ^  i  s  ^not  ^a^parameter  . 

end 

i  =  i  +  2 ; 

if  i  =  length  (  varargin ) 

error  (  ’No^name^  specified  ^for  „ parameter  .  ’  ) ; 

end 

end 


end 

function  r  =  Results  (obj) 
r  =  obj  .  Results  ; 

end 

function  disp(obj) 

fprintf(  ’Required:\n’  ); 
disp(obj  .  required  ); 

fprintf(  ’Optional;\n’  ); 
disp(obj  .  optional  ); 

fprintf  (  ’Named:\n  ’  ) ; 
disp  (  obj  .  params  ) ; 

end 


64 


UNCLASSIFIED 


UNCLASSIFIED 


DSTO-TR-2517 


B.6  Testing 

function  fail  =  t  es  t  _di  s  t  r  i  bu  t  i  on  ( id  ,  F,  samples,  w,  dist.alpha  ,  type) 

%  TEST.DISTRIBUTION  Test  the  Simulator’s  uncorrelated  distribution. 

% 

%  The  test-distribution  function  uses  Kolmogorov— Smirnov  to  test  the 
%  output  distribution. 

% 

%  Use  type  =  1  when  using  with  xUnit. 
ctx  =  Simulator  (F,  samples,  w) ; 

fail  =  0; 
i  f  type  =  0 

fprintf(  ’Test  u.%d^  (  uncorrelated  ^simulation  )\n’  ,  id); 
fprintf  (  ’ \tKolmogorov— Smirnov . ’  ) ; 

end 


X=  simulate  (  ctx  ) ; 
min_x  =  niin(X); 
max.x  =  niax(X); 


%  If  we  are  using  Matlab ,  then  we  use  the  builtin  kstest  funetion  . 
if  'exist  ( ’OCTAVE.VERSION’  ,  ’builtin’) 

cdf_points  =  min_x  :  (  (max_x— min_x  )  /  1 00  ):max_x; 


[ dist .fail  p] 


else 


kstest(X,  [cdf.points;  F( cdf. points )]’ , 
dist  .alpha  ) ; 


%  With  Octave  ,  things  are  somewhat  more  tricky  .  We  need  a  funetion 
%  entitled  <name>-cdf.  As  the  distribution  can  change,  we  have 
%  made  such  a  function  that  calls  the  function  — handle  in  a  global 
%  variable  id-cdf-f  ,  allowing  the  builtin  Kolmogorov— Smirnov  to 
%  work. 


%  If  id-cdf-f  is  already  set,  then  save  its  current  value  and  clear 
%  it  . 


UNCLASSIFIED 


65 


DSTO-TR-2517 


UNCLASSIFIED 


if  exist  (’ id_cdf_f  ’  ,  ’variable’) 
old-id.cdf. exists  =  true; 
old_id_cdf_value  =  id_cdf_f; 

i f  “isempty (whos(  ’global  ’  ,  ’variable  ’)) 
old-cdf -global  =  true; 

else 

old -cdf -global  =  false ; 

end 

clear  id  -cdf-f  ; 

else 

old-id-cdf-exists  =  false; 

end 

%  Set  our  CDF  and  run  K—S  on  the  dataset . 
global  i  d  -C  d  f -f  ; 
id-cdf-f  =  F ; 

p  =  kolmogorov-smirnov-test  (X,  ’identity’) 
if  p  <  dist-alpha 

dist -fail  =  true  ; 

else 

dist -fail  =  false  ; 

end 

%  Restore  id-cdf-f  to  its  old  value. 
if  old-id-cdf-exists 
clear  id-cdf-f  ; 

if  old -cdf-global 

global  id-cdf-f  ; 

end 

id-cdf-f  =  old -id -cdf-value  ; 
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end 

end 

%  Print  output  or  provide  success/ failure  status  to  xUnit. 
if  type  =  0 

if  dist -fail  =  I 

fail  =  fail  +  I; 

fprintf(  ’failed . ^ (p^=„%f ) \n\n  ’  ,  p ) ; 

end 

if  fail  =  0 ; 

fprintf(  ’passed. \n\n’  ); 

end 

else 

message  =  sprintf(  ... 

’  Kolmogorov— Smirnov^  failed  ^at^level  ^  alpha  ^=^%f  ’  ,  .  .  . 
dist -alpha  ) ; 

assertEqual  (  dist -f  ail  ,  false  ,  message); 

end 

end 

function  fail  =  t  es  t -Cor  r  e  1  at  i  on  ( id  ,  F,  r,  samples,  realisations, 

maxlag ,  w,  ... 
cdf-alpha  ,  ... 

cdf-error-threshold  ,  ... 

covariance-error-threshold  ,  ... 

test-type) 

%  TESTJJORRELATION  Test  the  Simulator  with  autocorrelated  data. 

% 

%  Use  test-typ  e  =  1  when  testing  with  xUnit. 

ctx  =  Simulator  (@(  t )  F(t),  r  ( 0  :(  samples  —  1 )) ,  w) ; 

if  test-type  =  0 
fail  =  0; 
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fprintf(  ’Test  (correlated^simulation)\n’,  id); 

end 

X=  simulate  (  ctx  ) ; 

if  exist  ( ’OCTAVE.VERSION’ ,  ’builtin’) 
cdf_points  =  unique(X); 

cdf-values  =  empiricaLcdf  (  cdf.points  ,  X); 

else 

output  =  zeros(l,  samples*  realisations  ) ; 
output  ( I :  samples  )  =  X; 

end 

%  We  need  to  store  the  output  covariances  (so  that  we  can  test  for 
%  autocovariance  error)  and  the  output  itself  (so  that  we  can  run 
%  the  chi-square  test  for  the  distribution. 
cov_out  =  zeros  (  size(  —  maxlag  :  maxlag  )) ; 

for  i  =  1 :  realisations 
X  =  simulate  (  ctx  ) ; 

if  exist  (  ’ OCTAVE ATIRSION ’ ) 

cdf-values  =  cdf.values  +  empir ical_cdf  (  cdf_points  ,  X); 
else  %Matlab 

%  Save  our  current  realisation  for  later. 
output  (  (  samples  *  i  + 1 ):(  samples  *(  i +1))  )  =  X; 


end 


cov.out  =  cov.out  +  xcov(X,  maxlag,  ’coeff’); 

end 

lags  =  —maxlag  :  maxlag  ; 
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cov_out  =  cov_out  /  realisations; 

%  If  we  are  using  Octave,  we  look  at  the  difference  between  actual 
%  and  empirical  CDF,  as  there  is  no  Chi— square  goodness  of  fit  built 
%  in  . 

if  exist  ( ’OCTAVE.VERSION’ ,  ’builtin’) 

max_cdf_error  =  max(  abs  (  cdf  .values  /  r  e  ali  s  at  ion  s  —  F(  cdf.points  ) ) ) ; 
if  max.cdf .error  >  cdf. error. threshold 
cdf.fail  =  true  ; 

cdf. error  .text  =  spr  intf  (’„  (max„  error  „%f )’  ,  max.cdf.err  or  ) ; 

else 

cdf.fail  =  false  ; 

end 


else 

%  Matlab  has  a  Chi— Square  goodness  — of— fit  test  built  in,  so  we  use 
%  this  instead  . 

[cdf.fail  p]  =  chi2gof ( output ,  ’cdf’,  F,  ... 

’ nbins ’ ,  10 ,  ... 

’alpha’,  cdf.alpha); 

cdf. error. text  =  sprintf  (  ’  ^  (p^=^%f )  ’  ,  p); 
cdf.fail  =  logical  (  cdf  .fail  ) ; 


end 

%  We  test  the  error  by  simply  comparing  the  autocovariance  co  efficients 

%  output  with  the  target. 

ac. error  =  abs(cov.out  —  r ( abs ( lags ) ) ) ; 

max. ac. error  =  niax(  ac. error  ) ; 

%  If  we  are  not  using  xUnit ,  print  the  test  output  ourselves. 
if  test. type  =  0 
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fprintf (  ’\tCDF„error . ’); 

if  cdf -fail 

fail  =  fail  +  1; 

fprintf ( ’  failed  .%s\n  ’  ,  cdf_error_text  ) ; 

else 

fprintf(  ’passed. \n’  ); 

end 


fprintf(  ’\tAutocovariance^error  ...  ’); 
if  max_ac -error  >  covariance-error-threshold 

fprintf(  ’failed  .„Greatest„error  „%f  .  \  n\n  ’  ,  ... 

max-ac-error  ,  covariance-error-threshold); 

fail  =  fail  +  1; 


return ; 

else 

fprintf(  ’passed. \n\n’  ); 

end 

else 

%  Otherwise  call  xUnit  functions  . 


assertFalse(cdf-fail  ,  ’Distribution^test^failed.’); 
assertTrue  (  max-ac-error  <  covariance-error-threshold  , 
’Autocovariance^error^too^large  .  ’) 


end 


end 
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Appendix  C  Tables  of  Numerical  Results 


We  have  made  numerical  estimates  for  the  coefficients  of  (26)  for  several  distributions. 


Exp(l) 

K(c  =  1,  =  1) 

n 

l(V',en)P 

rn{l) 

w(l) 

0 

1.000000 

1.000000 

2.467396 

1.532604 

1 

0.815765 

0.184235 

1.339585 

0.193019 

2 

0.177391 

0.006844 

0.188737 

0.004282 

3 

0.006685 

0.000159 

0.004139 

0.000143 

4 

0.000134 

0.000025 

0.000085 

0.000058 

5 

0.000017 

0.000008 

0.000004 

0.000054 

6 

0.000007 

0.000001 

0.000001 

0.000054 

Table  Cl:  Estimates  of  \{'ip,en)\‘^  for  the  Exponential  and  K(l,  1)  distributions. 


K(c  =  1 

,v  =  2) 

K(c  =  l 

v  =  A) 

n 

l(V’,en)P 

rn{l) 

w(i) 

0 

5.551652 

2.448348 

11.806900 

4.193100 

1 

2.235475 

0.212873 

3.939934 

0.253167 

2 

0.210615 

0.002258 

0.252644 

0.000522 

3 

0.002197 

0.000060 

0.000485 

0.000037 

4 

0.000054 

0.000006 

0.000024 

0.000013 

5 

0.000000 

0.000006 

0.000011 

0.000002 

6 

0.000000 

0.000006 

0.000002 

0.000001 

Table  C2:  Estimates  of  \{'ip^en)\‘^  for  the  2)  an(iK(l,4)  distributions. 


UNCLASSIFIED 


71 


DSTO-TR-2517 


UNCLASSIFIED 


N(^  =  I 

,(7  =  2) 

LN(//  =  0,  cr  =  1) 

n 

\{'4’,en)\^ 

rnil) 

rn{l) 

0 

1.000000 

4.000000 

2.71828 

4.67077 

I 

4.000000 

0.000000 

2.71828 

1.95249 

2 

0.000000 

0.000000 

1.35914 

0.59335 

3 

0.000000 

0.000000 

0.45305 

0.14030 

4 

0.000000 

0.000000 

0.11326 

0.02704 

5 

0.000000 

0.000000 

0.02265 

0.00439 

6 

0.000000 

0.000000 

0.00378 

0.00062 

Table  C3:  Estimates  of  \{^p,en)\‘^  for  the  Normal  and  Log-Normal  distributions. 
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